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We present a measurement of the atmospheric iZe, spectrum at energies between 0.1 TeV and 
100 TeV using data from the first year of the complete IceCube detector. Atmospheric Ue, originate 
mainly from the decays of kaons produced in cosmic-ray air showers. This analysis selects 1078 fully 
contained events in 332 days of livetime, then identifies those consistent with particle showers. A 
likelihood analysis with improved event selection extends our previous measurement of the conven¬ 
tional Ue fluxes to higher energies. The data constrain the conventional Ue flux to be 1.3lg'3 times a 
baseline prediction from a Honda’s calculation, including the knee of the cosmic-ray spectrum. A fit 
to the kaon contribution (f) to the neutrino flux finds a kaon component that is ^ = 1.3lg'4 times 
the baseline value. The fitted/measured prompt neutrino flux from charmed hadron decays strongly 
depends on the assumed astrophysical flux and shape. If the astrophysical component follows a 
power law, the result for the prompt flux is O.OjiQ'g times a calculated flux based on the work by 
Enberg, Reno and Sarcevic. 

PACS numbers: 95.55.Vj, 14.60.Lm, 29.40.Ka, 95.85.Ry, 25.30.Pt 

I. INTRODUCTION trino oscillation physics. Atmospheric muon and electron 

A measurement of the atmospheric neutrino flux is 
valuable in the field of neutrino astronomy and neu- 
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neutrinos are the decay products of mesons and muons 
which are produced when cosmic-ray primaries interact 
in the atmosphere. Experiments have measured the at¬ 
mospheric neutrino fluxes [ll-[l^ , and multiple theoretical 
frameworks to calculate this flux are available lll-ll7|. 

Below the knee (3 x 10^® eV) of the cosmic-ray energy 
spectrum, the flux of and Ve from tt and K decays, 
called the ‘conventional’ neutrino flux, follows a power 
law dN/dE oc where N and E are the number 

of neutrinos and the neutrino energy, respectively. The 
spectral slope is steeper than that of the primary cosmic 
rays by about one power because the neutrinos’ parent 
mesons lose a significant amount of energy in flight before 
decaying. 

The flux of conventional has been measured in a 
wide energy range. At energies below several lO’s of 
GeV, the flux is measured using fully contained events 
while, at energies above 100 GeV, flux measurements 
use muons produced by neutrinos traveling through the 
Earth, i.e. the upward-going direction. 

Most Ve come from the semileptonic decay of charged 
and neutral kaons. The Ve flux is lower than that of 
and the v^jv^. ratio increases with increasing energy, 
reaching a factor of ^20 at 1 TeV. The conventional 
and Ve flux is highest around the horizon, where parent 
mesons spend a higher fraction of their lifetime at higher 
altitudes and are less likely to interact before they can 
decay. 

The flux of high-energy conventional neutrinos is sen¬ 
sitive to the details of particle production in air showers. 
Large uncertainties on the conventional flux models at 
neutrino energy above 1 TeV come from uncertainties 
in strange quark production and the cosmic ray spec¬ 
trum, which are poorly constrained by accelerator and 
air-shower measurements. Precise measurements of the 
conventional and fluxes probe pion and kaon pro¬ 
duction in air showers. 

At energies between 1 TeV and 100 TeV, another class 
of atmospheric neutrinos arises, from charmed hadron 
decays. Since these hadrons have short lifetimes, the 
‘prompt’ neutrino flux retains the original spectral slope 
of the primary cosmic-rays. Prompt neutrinos are uni¬ 
formly produced in the atmosphere, with equal fluxes of 
and v^. The transition from the region dominated 
by the conventional neutrinos to the prompt neutrinos in 
the spectrum is expected to occur at energies of around 
1 PeV for and around 30 TeV for v^. 

Theoretical predictions for the atmospheric charm pro¬ 
duction have large uncertainties l3, largely due to 
a lack of data on forward production at high energies. 
RHIC and LHC data provide useful constraints, but only 
in the central region. Several non-perturbative effects 
come into play in the forward region of collisions. Uncer¬ 
tainties in the low—a: parton distributions and possible 
diffractive production channels lead to significant uncer- 


mrr* 

tainties |19^21|. 


Another flux component is the astrophysical neutri¬ 
nos recently discovered by IceCube [2^ |23[ . At energies 
above lO’s of TeVs, a seemingly isotropic flux of neu¬ 
trinos of astrophysical origin becomes discernible with 
a spectrum harder than that of the atmospheric flux. 
However, it is difficult to disentangle the prompt flux 
from the astrophysical component with the current event 
samples because their angular distributions and spectral 
indices are similar. Recent IceCube analyses address the 
issue 


24, 


In this paper, we present a measurement of the atmo¬ 
spheric r'e spectrum with IceCube. 


II. DETECTOR 


IceCube is a neutrino detector deep in the South Pole 
ice. The cubic-kilometer detector consists of 5,160 light 
sensors distributed on 86 vertical strings at depths be¬ 
tween 1450 m and 2450 m below the surface. The array of 
sensors, called Digital Optical Modules (DOM) [i^, ob¬ 
serves Cherenkov light produced when charged particles 
exceed the speed of light in the ice. The DOMs consist 
of a pressurized glass sphere, a 252 mm diameter Pho- 
toMultiplier Tube (PMT) and digitizing electronics. 
Twelve LEDs in each DOM are used to calibrate the de¬ 
tector responses. 

The calibration of the DOM response and the under¬ 
standing of the optical properties of the surrounding ice 
are crucial for the event reconstruction in IceCube. Using 
in situ LED data, the ice is modelled as a set of scatter¬ 
ing and absorption parameters as functions of wavelength 
and depth [28|, l29j. The ice exhibits an optically layered 
structure depending on dust concentration, reflecting the 
long-term differences in climate that affected dust accu¬ 
mulation over time. 

The IceCube neutrino observatory includes three com¬ 
ponents, each designed for a specific purpose. The 
baseline array contains 4,680 DOMs on 78 strings with 
roughly 125 m string-to-string distance and 17 m DOM- 
to-DOM spacing and is optimized for detecting neutri¬ 
nos above a few 100 GeV. The “DeepCore” sub-array 
is a more densely instrumented set of DOMs optimized 
for identifying neutrino events with energies as low as 
10 GeV [^. It contains 480 DOMs on 8 strings de¬ 
ployed in the bottom-center part of the baseline array 
together with DOMs of the baseline array in the same re¬ 
gion. Air showers are observed by a surface array called 
IceTop [m. 

The DOMs digitize the recorded PMT waveforms and 
generate time-stamped signals, or “hits” when the sig¬ 
nal rises above a threshold which is set to 0.25 pho- 
toelectrons. The ATWD (Analog Transient Waveform 
Digitizer [s^) and fADC (fast Analog to Digital Con¬ 
verter) digitize the waveforms at rate of 300 and 25 mega¬ 
samples/s, respectively. The ATWD records 128 samples 
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(430 ns total) with a charge resolution of ~30% for sin¬ 
gle photoelectrons and a timing resolution of ^2 ns. The 
fADC system records 256 samples/event (6400 ns), to 
capture long, late pulses. If a nearest or next-to-nearest 
neighbor DOM is also hit within ±1,000 ns, then the 
DOM transmits the full waveforms to the surface. Oth¬ 
erwise, for isolated hits, it sends a brief summary. The 
isolated hits are important for efficiently recognizing in¬ 
cident cosmic-ray muons which might give a faint light 
through minimum ionization. 

The surface electronics forms a trigger when at least 
eight non-isolated hits are observed in a 5,000 ns win¬ 
dow. Then, a physics event is built, containing all of the 
isolated and non-isolated hits. Further details about the 
detector can be found in Ref. [ii-Isi. 

The IceCube coordinate system is right-handed, with 
its origin at the center of the baseline array, with the 
z-axis pointing upward. The y-axis follows the prime 
(Greenwich) meridian, and the x-axis points toward ±90 
degrees longitude. The zenith angle (9) is defined in the 
usual manner, the angle between the event arrival di¬ 
rection and the z-axis, while the azimuthal angle ((/>) is 
measured from the positive cc-axis, in the x — y plane. 


III. DATA AND SIMULATION 

This analysis uses data taken with the full 86-string 
configuration of IceCube, between May 13, 2011 and May 
15, 2012. After excluding calibration runs and a few pe¬ 
riods when the detector was operating in a partial config¬ 
uration or exhibiting large variations in rate, the livetime 
is 332.3 days. This selection avoids systematic biases due 
to detector instability and ensures all strings of the de¬ 
tector are active. 

In order to avoid statistical bias, the analysis cuts and 
fit procedure were developed using only 10% of the data, 
spread evenly throughout the year. After the cuts and 
fit were fixed, the rest of the data was studied. 

For this analysis, the signal is defined as atmospheric 
interactions contained inside the detector volume. Con¬ 
tainment criteria are based on the vertex position, which 
is determined using both the first DOM hit in time and a 
vertex reconstruction. Non-contained background events 
entering from outside the detector are vetoed by these 
containment criteria (Section IV). 

When a high-energy neutrino interacts in the ice, the 
deep inelastic scattering can produce one of three event 
signatures. ‘Cascades’ are created by Ve charged cur¬ 
rent (CC) interactions which consist of an electromag¬ 
netic shower and a hadronic shower, or neutral cur¬ 
rent (NC) interactions of all neutrino flavors. ‘Tracks’ 
are through-going muons from CC interactions occur¬ 
ring outside of the detector. ‘Hybrid’ events from CC 
interactions occurring within the detector have both a 
hadronic shower and a track. 


The Cherenkov light yield of the shower particles is 
proportional to the cascade energy. Hadronic showers 
have lower light output and larger shower-to-shower vari¬ 
ations than electromagnetic showers 36, [s^. This is 
partly because hadrons are heavier than electrons, with 
higher Cherenkov thresholds. Also, hadronic showers 
produce neutral particles, have nuclear interactions, and 
transfer energy to struck nucleons. A 1 TeV hadronic 
shower has a light output which is (80± 10)% [s^ of that 
of an electromagnetic shower of the same energy . In 
the simulations, a parameterization is used to account for 
the reduced light output. The visible energy (Emg) is de¬ 
fined as the observed energy, assuming that the shower is 
electromagnetic. The pattern of detected light is roughly 
spherical for both types of showers due to short travel 
lengths of the shower particles. 

The largest background in this analysis is downward 
going muons produced by high-energy cosmic ray inter¬ 
actions (‘CR muons’). The CR muons which reach the 
surface with an energy of 500 GeV or larger can penetrate 
the ice to the depth of IceCube, and can become a back¬ 
ground to the atmospheric neutrino signal. This muon 
background has three main signatures. The first kind is 
from through-going tracks created outside of the detec¬ 
tor. This could be a down-going single muon or muon 
bundle from a cosmic-ray interaction. The second kind 
is an event with multiple tracks having different direc¬ 
tions produced by coincident but unrelated air showers. 
The third kind is a ‘stealth’ muon which passes between 
strings or through the dustiest, optically most absorb¬ 
ing ice layers. This class of events has the appearance of 
the cascade signal when the muon generates only a few 
hits in the outer regions of the detector and then under¬ 
goes stochastic losses that release most of its energy in a 
cascade-like shower within the fiducial region. There is 
also a small background from through-going muons from 
interactions outside of the detector. 

Air showers are simulated with CORSIKA (COs- 
mic Ray Simulations for KAscade) [39| including the 
Sibyll |40| hadronic interaction model. In IceCube the 
cosmic ray spectra are simulated for five nuclei. By re¬ 
weighting the five spectra, a resulting muon flux is ob¬ 
tained to represent a cosmic ray composition model. In 
this analysis, we used the phenomenological ‘H3a’ com¬ 
position model ^ which takes into account updated cos¬ 
mic ray spectra and the most recent spectral slope mea¬ 
surements |42l444l| . IceCube data are in good agreement 


with the H3a model predictions in the energy range rele¬ 
vant for this analysis (1 — 1000 TeV in the primary cosmic- 
ray energy). An alternative model, the poly-gonato spec¬ 
trum [45l| also models cosmic rays with five different nu¬ 
clei. It uses different parameterizations, particularly for 
the knee, and finds, for this analysis, a roughly 30% dif¬ 
ference in the background normalization. 

Neutrinos are simulated by using a software mod¬ 
ule (neutrino-generator) which is based on the ANIS 








5 


package [46| with CTEQ5 cross-section tables. To 
obtain a large sample of simulated neutrinos, the simu¬ 
lation forces all neutrinos to interact near the detector 
and then assigns each an interaction probability. The 
conventional flux from Ref. ll| (‘Honda’) was used up 
to 10 TeV. The Honda model is extrapolated to higher 
energies using the flux parameterization 


‘^{E,) = C -E-^ -{w^+wk), 

A 


W-jr = 


I + COS 9*/eTr’ 

_ A-Kv _ 

1 -I- BkvEv cos 6*!ex' 


WK = 


( 1 ) 

( 2 ) 

(3) 


The Wjr and wk are relative contributions to the neutrino 
flux from TT and respectively. The parameters H, B, 
and the absolute normalization (C) of the flux are de¬ 
termined by fitting to the published Honda flux at lower 
energies and 9* is the neutrino zenith angle at the produc¬ 
tion point. The index and critical energies are a = 2.65, 
e,r = 115 GeV, and ey = 850 GeV. For prompt neutrinos, 
the flux from Ref. [l^ (‘ERS’) is used. Both the conven¬ 
tional and prompt baseline predictions are corrected to 
an updated cosmic-ray spectrum, including a knee struc¬ 
ture which is similar to the H3a spectrum used for the 
cosmic-ray simulations. We also apply a small correction 
factor of 0.5% to account for the additional I'e produc¬ 
tion from Ks semileptonic decays (d^ . Tau neutrinos are 
not included as an atmospheric component. The con¬ 
tribution is less than 5% compared to the total prompt 
contribution because its flux mainly comes Dg decays 
which are smaller than other charmed hadron decays in 
the atmosphere M- 

With increasing energy, the probability of vetoing an 
atmospheric neutrino through the presence of CR muons 
from the same cosmic-ray shower increases in the down¬ 
ward region. Since the two coincident particles are nearly 
collinear, the events are automatically rejected in anal¬ 
yses sensitive to the downward contained events and 
therefore the veto probability as an additional correction 
should be applied to the event rate 0,[5ll. 

The efficiency of this atmospheric self veto correction 
depends on the neutrino production processes. There 
are two types of veto. A ‘correlated’ veto occurs when 
a muon and a neutrino are produced in the same decay, 
while for a ‘non-correlated’ veto they have different par¬ 
ent particles from the same shower. Eor the conventional 
the correlated veto is the dominant process while, 
for the conventional iZe, the non-correlated veto is the 
main process. Eor the prompt neutrinos at energies well 
above the ex, the impact from the non-correlated compo¬ 
nent increases since decay modes of the parent charmed 
hadrons involve fewer number of correlated muons. The 
reference model [sH used here treats the correlated and 
non-correlated components for and iZg separately. As 
Fig. 3 of Ref. ^ shows, for neutrinos above 10 TeV 


in the vertically downward region, the veto probability is 
higher than 95% while at cos{9) = 0.2, the probability is 
close to 50%. 

In the following, the ‘modified Honda’ flux includes the 
extrapolation to higher energies of the original Honda 
flux with the input H3a spectra and the additional Ks 
contribution. Similarly, the ‘modified ERS’ flux refers to 
the ERS flux with the H3a spectra. These modified fluxes 
are used to weight the neutrino simulations throughout 
this analysis (details can be found in [I^.) For base¬ 
line event rate predictions, the self veto corrections are 
applied to these modihed fluxes. To ensure consistency, 
the conventional neutrino expectations are partly vali¬ 
dated against the neutrino predictions from the COR- 
SIKA generator in the downward region where the self 
veto is in effect. The astrophysical flux is modelled as a 
single power law dN/dE ^ (j)o ■ (E/lOOTeV)”'’', where (j)o 
is the flux at 100 TeV and 7 is the spectral index. 

The simulated particles are propagated to the detector, 
and Gherenkov photons are produced from the charged 
particles. The generated photons are tracked through 
the ice, using the measured optical scattering and ab¬ 
sorption coefficients, and then through a simulation of 
the hardware response. For example, at z = —350 m, 
a Cherenkov photon can travel 200 m with on average 
three scatters. IceGube simulates photons from an event 
in a cylindrical volume which is roughly twice bigger in 
radius and length than the instrumented volume to max¬ 
imize the light collection efficiency for events skimming 
the edge of the detector. The Monte Carlo events use the 
same format as data events. Both are treated identically 
in event processing. The average trigger rate was 2200 
Hz with a roughly 10% seasonal variation due to tem¬ 
perature and pressure changes in the atmospheric condi¬ 
tions above IceCube. The CR muon rate at trigger level 
(‘Level 1’) is 7.3 x 10® times larger than that of the at¬ 
mospheric i/g signal, which is predicted to be 3.0 x 10“^ 
Hz above 300 GeV. 


IV. EVENT SELECTION 


The event selection proceeds in several stages to enrich 
the atmospheric Ve signal against the large GR muon 
background by contrasting the simulated signal with 
background Monte Carlos. As it is currently not possi¬ 
ble to distinguish electromagnetic showers from hadronic 
ones in IceCube, the sample contains background from 
NC interactions of other neutrino flavors. The selection 
relies on the searches for the spherical hit pattern of light 
in the cascade signal and reconstruction variables which 
describe the cascade signal in ice. The cascade variables 
used in this analysis are explained in more details in 
Ref. dm. 
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Reconstructions 

Two maximum-likelihood algorithms are used to re¬ 
construct events under the cascade hypothesis. The Hrst 
estimate (‘Cascade-LLH ’) [s^, of a cascade inter¬ 
action position (.Aivertex: k^ertexi and .^vertex); and an 
interaction time uses only hit time information. The 
algorithm uses an analytic probability density function 
(PDF), p{tres,di) [ 5 ^ expressed in inverse nanoseconds, 
for constructing a likelihood : 

hits 

C = Y\_P{tres,di) (4) 

tres — igeo — tg di / Cice^ (5) 

where is the observed time of the hit, tg the ex¬ 
pected time of the cascade interaction, di the distance 
from the hit DOM to the interaction vertex, and Cice the 
speed of light in ice. The time delay of a hit relative 
to the geometrical time (tgeo) corresponding to straight- 
line propagation is dehned as the residual time, t^es, he. 
a non-scattered photon registers at t^es = 0. Cascade- 
LLH provides an initial vertex seed for an improved re¬ 
construction and returns a cascade quality parameter 
RLLHyertex = —^ogC/{Nhit — 4), an analog of reduced 
X^, with a number of hits {Nhu) minus four degrees of 
freedom. The calculation assumes that photon scattering 
is independent of depth in the ice, so the reconstruction 
has a limited resolution. The resolution is measured as a 
Gaussian spread of 1 tr on the difference between a true 
cascade vertex and a reconstructed vertex. The Cascade- 
LLH reconstructs the interaction vertex with a resolution 
of 11 m in the x-y plane and 12 m in z for 10 TeV v^. 

The second, more advanced algorithm, (‘CREDO ’) [l|, 
[sol . [ 5 ^ reconstructs seven parameters of a cascade in a 
single fit. The vertex position (X^eco,Yreco, and ^reco), 
the time (treco), the direction {0reco and 4>reco) and the 
visible energy {Ereco) of the cascade are estimated by 
using full waveform information. CREDO uses a more 
detailed PDF for time and amplitude expectations that 
includes the depth dependent propagation of light in ice. 
The scattering and absorption properties are stored in a 
table which is interpolated with splines [^. The vertex 
resolution of CREDO is 4 m in the x-y plane and 3 m in 
z for 10 TeV i/g. 

Likelihood reconstructions based on the track hypoth¬ 
esis [H^ are used to identify the CR muon background 
events and to estimate their direction {0track and (fitrack)- 
As a muon travels close to the speed of light, a likelihood 
similar to Eq. |T]with a new tgeo definition is fitted. The 
PDE used in the reconstruction is called a single pho¬ 
toelectron (SPE) PDE which models the tres using only 
the earliest photon at each DOM. Additionally, a zenith- 
weighted Bayesian track reconstruction is performed us¬ 
ing prior knowledge of the CR muon angular distribu¬ 
tion. Only the downward-going direction is allowed for 


the reconstructed track directions since the reconstruc¬ 
tion maximizes the product of the PDF and the prior. 

Level 2 (Cascade Online Filter) 

The events recorded at the South Pole are filtered to 
reduce the data volume so that the data can be trans¬ 
ferred to off-sites via satellite. This online filter (Level 2) 
algorithm removes early and late PMT hits unrelated to 
physics interactions. The remaining hits are used to cal¬ 
culate simple topology variables and RLLHyertex which 
are used in the filtering. The background rejection factor 
is 99% with the filter retaining 77% of the atmospheric Ve 
signal above 300 GeV. The efficiency reaches 90% above 
10 TeV. The CR muon background after Level 2 selec¬ 
tion comprises about 60%, 20%, and 20% of through- 
going muons, coincident muons, and stealth muons, re¬ 
spectively. 

Level 3 (Containment) 

The containment cuts require that the cascade vertex 
is in the fiducial region. In addition to simple contain¬ 
ment conditions based on the earliest hit time, cuts based 
on the vertex reconstruction are applied to make the light 
produced by the cascades contained within the detector 
volume. This Level 3 filter reduces the CR muon back¬ 
ground further as the background-to-signal ratio is still 
high (~ 10'^). 

An algorithm identifies clusters of hits which are dis¬ 
tinct in time and space. Only the events classified as a 
single cluster are accepted, in order to reject coincident 
CR muon background. 

The first hit must not be on one of the outer strings, 
and must be no closer than 70 m to the top or bottom of 
the detector, i.e. —430 m < Zi < 430 m. 

The fiducial volume cut requires the reconstructed ver¬ 
tex from CREDO (Vreco, V-eco, Zreco) must be within a 
cylinder of 420 m radius from the center of the detector 
with 70 m minimum distance to the edge of the detec¬ 
tor. Additionally, the vertex should be no closer than 
100 m to the top or 50 m to the bottom of the detector 
(—450 m < Zreco < 400 m) . 

For the contained events, we further impose several 
quality cuts. Each event must have hits on at least three 
non-DeepCore strings. The ratio of the number of hit 
DOMs to the total number of DOMs within a ^here 
centered on the vertex should be high (>60%) [l|, l59l| . 
The radius of the sphere is determined by the root-mean- 
square distance to the vertex of hit DOMs with a scal¬ 
ing factor that maximizes signal selection power. Then, 
events with a low RLLHyertex are selected. After this 
Level 3 selection, the dominant CR muon background 
is stealth muons with a few veto hits. These muons 
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are typically minimally ionizing in the veto region and 
then produce a stochastic signature in the fiducial re¬ 
gion, mimicking the cascade signal. The cut efficiencies 
will be discussed later. 


Level 4 (Neutrino Selection) 


TABLE I. The number of events in 332 days of data are shown 
at each level. The Level 4 rates are the final best fit values 
where CR muons are estimated from the data side band. The 
total neutrino rates are based on the modified Honda predic¬ 
tions. Level 1 & 2 neutrino rates are only for energies above 
300 GeV. i^e contains about 8% NC events at all levels but 
NC fraction(J?jvc) grows as the higher level selection enriches 
cascades. 


The Level 4 event selection uses a machine learning 
technique to separate the atmospheric cascade signal 
from the CR muon background. A multivariate analysis 
method based on boosted decision trees (BDT) is imple¬ 
mented using a toolkit for multivariate data analysis [60j . 
The BDT uses 12 variables, chosen for their power to 
separate cascades from the CR muon backgrounds. The 
variables are listed below, classified in three categories: 
veto, quality, and topology. 

Veto: 

1- Nyeto'- The number of hits recorded before 
the CREDO vertex in time and consistent with 
downward-going muons. 

2. Ncone- The number of hits in a cone with its 
apex at the CREDO vertex position, an opening 
angle of 36 degrees, and centered on the incoming 
track direction obtained from the SPE reconstruc¬ 
tion. 


Preco \/^reco T ^eco 


Quality: 

5. 9track'- The zenith angle of the SPE reconstruc¬ 
tion. 


6 - dreco- The zenith angle of the CREDO recon¬ 
struction. 


7. RLLHy 


8 . Riikeiihood- The likelihood ratio of the SPE- 
track hypothesis to the Cascade-LLH hypothesis. 
9- Rsayes- The likelihood ratio of downward 
forced track hypothesis (zenith-weighted SPE) to 
non-forced track hypothesis (SPE). 


Topology: 

10. Rq: The charge fraction in the first 300 ns, 
excluding the two earliest hits. 

11. Zsplit- The vertical distance between hit 
center-of-gravities determined by splitting in time 
into two clusters of hits. 

12. Zspeed- The z-coordinate component of a 
reconstructed velocity calculated using the first 
half of all hits in an event. 


The discriminating power comes relatively evenly from 
the three categories. The most powerful separators are 


Level 

Data 

CR muon 

Atm. Ve[RNC%) 

Atm. Vp,{Rnc%) 

1 

6.3 X 10^° 

5.2 X 10^° 

8.6 X 10® (10%) 

2.3 X 10® (6%) 

2 

7.8 X 10® 

5.8 X 10® 

6.6 X 10® (8%) 

1.0 X 10® (14%) 

3 

2.3 X 10® 

2.3 X 10® 

1278(8%) 

5784(35%) 

4 

1078 

115 

215 (8%) 

645 (40%) 


Rukeiihood, zsplit, and Ncone which not only describe the 
data well but also show minimal correlations with other 
variables. The distributions of these variables are shown 
in Fig. [TJ The Veto variables ensure that the contained 
events show no trace of an incoming muon track before 
a reconstructed vertex time. The Quality variables iden¬ 
tify whether the hit pattern is extended in any particular 
direction (track-like) or is isotropic (cascade-like). The 
Topology variables look for cascade events with more lo¬ 
calized charge distribution to distinguish them from the 
long through-going events that distribute hits in a larger 
distance. The BDT output is a discriminant score for 
each event, with a higher number indicating more signal¬ 
like events. 

The distributions of the BDT scores are shown in 
Fig. [21 The data and the total Monte Carlo prediction 
display a transition from the region dominated by the 
CR muons to the atmospheric neutrino dominant region. 
The gradual over-prediction of the CR muon background 
going from a low BDT score to a high score is mainly 
due to the limitations in the modelling of the detector 
systematic uncertainties. At high BDT scores, the CR 
muon background events are relatively more populous at 
the bottom part of the detector where our veto is less ef¬ 
fective in rejecting the CR muon simulation events than 
in data. These events come through the dustiest ice re¬ 
gion between 2000 m and 2100 m in depth with shallow 
zenith angles and produce few hits in the veto region. 
Since the BDT uses depth-dependent variables such as 
Zreco, it is sensitive to absorption of the photons in ice. 

The final selection is based on the two-dimensional 
BDT-Energy cut shown in Fig. |3l One limitation of this 
selection is that it was based on CR muon simulations 
which had limited statistics (about 10% of the data live 
time), set by the available processing power. For this rea¬ 
son, the CR muon background was not estimated from 
the CORSIKA sample using the final BDT-Energy cut. 
Table U summarizes the event selection. The efficiency is 
shown in Fig. U) 
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FIG. 1. Four representative variables in the BDT (see description in the text). The dotted red lines are for CR muon background 
simulation while the blue lines are for the atmospheric Ve. simulations. Real data are shown with black circles. The signal is 
scaled up by 1000 to show the separation. 



FIG. 2. BDT score distribution at Level 3. Real data are 
shown with points while the sum of all Monte Carlo predic¬ 
tion is shown as a solid green line. Atmospheric neutrino 
predictions are shown with the cyan dotted line (i'^) and the 
blue solid line [ve)- The CR muon background simulation is 
shown with a red line. The side band for the final level CR 
muon estimation is indicated with the magenta vertical band. 


To get a more accurate estimate of the CR muon back¬ 
ground without relying on simulations, we estimate the 
final rate from background-dominated data close to the 
signal region. As Fig. [3] shows, the region around 10 TeV 
is poorly populated in the CR muon background sim¬ 
ulation. However, the BDT score shows no strong de¬ 
pendence on the reconstructed energy. This is expected 
because the training variables do not contain explicit en¬ 
ergy information. The data region for the background 
estimation is chosen such that sample size is maximized 
and neutrino contamination is minimized while staying 
as close as possible to the high-BDT score signal region. 
We obtain the optimal control sample with BDT scores 
between 0.25 and 0.3, as shown by the vertical band in 
Fig. [3] Observable distributions for data in this region 
are used as templates for the final CR muon background. 
To check the robustness of this background estimate, 
distributions of neighboring data regions with the same 
width have been evaluated. Results using these alterna¬ 
tive bands show no significant deviation from the baseline 
choice. 


V. PARTICLE IDENTIFICATION 

After Level 4, particle identification (PID) variables 
are used to distinguish hybrid (v^ CC) events from cas¬ 
cades. For a more accurate vertex position, angular, and 
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FIG. 3. Two-dimensional distributions of a reconstructed visible energy as a function of the BDT score are shown for the real 
data (left), CR muon background (middle), and atmospheric neutrinos (right). The right side of the black vertical line indicates 
the final event sample. The z-axis is the number of events in 332.3 days. 



log^^(EyGeV) 


FIG. 4. The Ve. effective volume is plotted as a function of 
energy at different cut levels. The effective volume is the de¬ 
tector volume multiplied by the ratio of the number of selected 
events to the number of generated events over 4 tt steradian. 
No self veto correction is applied to the number of events. The 
relative efficiency change at high energies is mainly due to con¬ 
tainment requirements while, at low energies, the change is 
driven by the strict selections needed to reject backgrounds. 


energy reconstruction of cascades, we use slower ‘iter¬ 
ative’ CREDO reconstruction P, [^, . The CREDO 

reconstruction with four different angular seeds mitigates 
the probability of the minimizer becoming trapped in lo¬ 
cal minima, a common problem in scenarios with a high 
number of dimensions. The improved results for a re¬ 
constructed energy {Ereco,A) and for the reconstructed 
zenith angle (0reco,4) are shown in Fig. [S] and used in the 
analysis fitting procedure. 

For signal CC events at around 10 TeV, the en¬ 
ergy resolution is AEyis/Eyis Ri ±9%, where AEyis = 
Eyis — Ereco,4- The mean of the AEyis distribution over¬ 
estimates the visible energy by about 6%, but, because 
of neutral-current interactions, Eyis is on average 4% 


lower than the true neutrino energy. For the same energy 
range, the zenith angle accuracy is AQyjQy k, ±8 degrees, 
where By is the true zenith angle of the neutrino and ABy 
is the difference between the neutrino zenith angle and 
Breco, 4 - The mean of the ABy distribution does not show 
a bias. 

Systematic effects add additional uncertainties on the 
resolutions which are evaluated using alternative simula¬ 
tions. The optical efficiency of a DOM and the optical 
properties of ice are varied for those simulations by a 
known amount. We treat a maximum deviation from the 
baseline simulation as a size of the uncertainty. From 
this study, 12% for energy uncertainty and 2 degrees for 
zenith angle uncertainty are obtained. 

With the iterative CREDO results, four selected vari¬ 
ables are used to train a BDT for particle identification 
(PID-BDT) by treating CC events as background. 
These variables exploit hits originating from the muon in 
the CC interaction. Two variables depend on the first 
photon arrival times at the DOMs, relative to the time 
expected for a point-like emitter (cascade hypothesis) for 
a photon that does not scatter in the ice. Since muons 
move faster than photons in the ice, they are likely to 
produce acausal photons. The residual ty^s is calculated 
for each photon. The variables are 

PID variables: 

Tl: The smallest ty&a 

T2: The number of hits in —200 ns < tres < 20 ns 
T3: The distance that the cascade vertex moves 
when it is reconstructed, after omitting the acausal 
hits from the reconstruction. 

X4: RLLEly^ytyx 

The PID-BDT output shown in Fig.[7]agrees well with 
the simulation expectation and shows good separation 
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FIG. 5. Cascade reconstruction uncertainties as a function of 
true Ue energy using Vs CC interactions. The upper panel 
shows the energy statistical (only) resolution obtained by 
comparing the visible energy with reconstructed energy. The 
lower panel is for zenith angle resolution in degrees. Solid 
blue lines represent the iterative reconstruction results. Sys¬ 
tematic uncertainties (gray bands) add an additional 12% for 
energy uncertainty and 2 degrees for zenith angle uncertainty. 



FIG. 7. PID-BDT output distribution for the baseline best 
fit. The conventional CC as a hybrid component and the 
conventional cascade component including i/g and CC are 
plotted separately. Events with high PID-BDT scores are 
cascades. The prompt component is fit to zero. 



FIG. 8. Reconstructed zenith angle for the baseline best fit. 
The prompt component is fit to zero. 


FIG. 6. Reconstructed energy distribution for the baseline 
best fit. Note that the prompt component is fit to zero. 

between CC and izg-like events. Monte Carlo studies 
show that the CC identification improves at higher 
energies as the muon track becomes more visible. 

VI. FLUX MEASUREMENT METHOD 

For measuring the atmospheric neutrino fluxes, the 
data was histogrammed in three dimensions — energy, 
zenith angle and PID observable. Three different fits 
were performed to the final data sample, to test differ¬ 
ent physics parameters. The first (baseline fit) was a 
straightforward measurement of the Vg spectrum, assum¬ 


ing the spectral shape of the model. The second subdi¬ 
vided the energy spectrum, to make a binned measure¬ 
ment of the r'e flux versus the neutrino energy. The third 
fit was similar to the baseline fit, but allowed the kaon 
to pion ratio to vary. These fits include the systematic 
uncertainties using a profile likelihood approach. 

The parameters of the baseline fit are shown in Ta¬ 
ble HU Six physics parameters are used: conventional 
and Vg normalizations relative to the modified Honda 
flux, a CR muon normalization, a total prompt (zz^ + Vg) 
normalization with respect to the modified ERS flux, 
an astrophysical normalization (t/fo) and an astrophysical 
spectral index ( 7 ). Figures 6 to 8 show one-dimensional 
projections of these histograms with the bin numbers and 
their ranges used in the fitter, along with the baseline fit 
results. 
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TABLE II. Baseline fit result table. The best fit includes 
statistical and systematic uncertainties at 68 % C.L. 


Parameters 

Best fit 

CR muon 

115^27 events 

Conventional 

l.Olg j X modified Honda 

Conventional Ve 

l.Slg’g X modified Honda 

Prompt Normalization 

O.Olg;^ X modified ERS 

Astrophysical flux 

3.2lJ;i X 10-1® 


xGeV-icm-%r-is-i(E,/ 10 ®GeV)-^ 

Astrophysical 7 


Optical Efficiency 

- 2 . 2 % (± 10 % prior) 

Ice Parameters 

± 1.0 a at scattering ± 10 % 


In the second fit, the conventional v^. flux component 
is further divided into four smaller energy ranges span¬ 
ning 100 GeV to 100 TeV, introducing three additional 
physics parameters. Because the region above 100 TeV 
is dominated by the astrophysical component, it is not 
used in this separate fit. 

Finally, in the third fit, a kaon fraction parameter and a 
total conventional [y^ -\-Ve) normalization are introduced 
in order to remove any correlation between conventional 
and Vf.. 

For these fits, the likelihood L is constructed with a 
Poissonian component for the physics parameters and 
Gaussian components for the systematic parameters. 
Best fit results are obtained by minimizing the negative 
logarithm of L, 


-21nL = 2 y^(^k-nkln/rk)(6 ) 

k m 

All physics parameters are unconstrained in the fitting 
process while the two systematic parameters are re¬ 
stricted by the priors that quantify their estimated preci¬ 
sion. The total expected count /x is a sum of each compo¬ 
nent’s contribution which depends on the fit parameters. 
The index k iterates over the histogram bins, and the 
number of observed events in bin k is rik. The expected 
count also depends on given systematic parameters. The 
two systematic parameters (m) have a central value {1) 
and an uncertainty (cr;). The 68% parameter uncertain¬ 
ties are determined by scanning — 21nL up to one unit 
from the best fit likelihood. 


VII. SYSTEMATIC UNCERTAINTIES 

Systematics uncertainties arise due to imperfect mod¬ 
elling of our detector which can affect analysis results. 
The two most important detector systematics are in¬ 
cluded in the fitter. They are the total optical efficiency 
of a DOM and the optical properties of surrounding ice 


■o 


- Nominal 

. DOM eff. +10% 

- DOM eff.-10% 


4 4.5 5 

log (E /GeV) 

ft ' vann & ' 



FIG. 9. Effects of systematic variations in final level Ve. sam¬ 
ple. The upper panel shows the DOM efficiency impact on 
the reconstructed energy compared to the nominal predic¬ 
tion. The lower panel shows the ice systematic impact due to 
an increased scattering on the reconstructed zenith angle. 


(scattering and absorption lengths). The sizes of these 
systematic errors are estimated from laboratory measure¬ 
ments of DOMs and ice measurements with in-situ de¬ 
vices. A few simulations with different input assumptions 
on the systematic effects are performed and their event 
rates and shapes are compared with those in the nominal 
Monte Garlo at the final analysis level. 

Simulations with modified optical DOM efficiency re¬ 
sult in a different event rate globally but show little 
change in the shape of the analysis observable distribu¬ 
tions. Since the normalization of an assumed physics 
model translates directly to a flux of that model, this 
systematic uncertainty loosens the constraint on the flux. 
The impact on the event rate relative to the nominal 
value is parameterized using five simulations with the in¬ 
put efficiencies ranging from —10% to -1-10%. This results 
in output event rate changes in the range [—20%, -|-10%], 
matching the high event rate for the high efficiency in¬ 
put (see Fig. HD. The asymmetric change in event rate 
with varying optical efficiency is a combined effect of 
the changing number of observed photoelectrons, coupled 
with the analysis selection of higher quality events. 

The ice systematics alter the shape of the zenith angle 
distribution. A global increase in light scattering tilts 
the zenith angle reconstruction downwards (see Fig. HD. 
This is a consequence of losing non-scattered hits that 
are crucial in the cascade direction reconstruction. A 
10 % increase in scattering coefficient degrades the angu¬ 
lar resolution by about 2 degrees. However, the change 
in the absorption coefficients has a smaller impact on the 
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zenith angle shape. The change in zenith angle distribu¬ 
tion is modelled as a parameterization of the ice model by 
reweighting the event rate. This model changes the nom¬ 
inal event rate by — 10 % in the upward going direction 
and by -|- 10 % in the downward going region for a positive 
one-sigma shift in the fitter (— = 1) (see Fig. |9l). 

We have investigated other systematic effects arising 
from the neutrino-nucleon cross-section and cosmic ray 
spectral slope. The theoretical uncertainties from the 
neutrino-nucleon deep inelastic cross-section ( 2 ^ [ 2 lj| are 
relatively small compared to the other systematic un¬ 
certainties. We assume a 3% cross-section uncertainty, 
following Ref. [ij. The systematic impact of the cross- 
section acts as a simple normalization in the energy re¬ 
gion of this analysis and is strongly correlated with the 
DOM efficiency parameterization. The cosmic ray spec¬ 
tral slope has a small impact, compared to the detector- 
related systematic uncertainties. Additionally, the sys¬ 
tematic effect on the zenith angle shape due to a seasonal 
temperature variation and the atmospheric self veto cal¬ 
culation are similar to ice systematics and absorbed by 
the parameterization of ice systematics. 


VIII. RESULTS AND DISCUSSION 

A total of 1078 events are observed after unblinding 
the full data set. The cascade candidates are distributed 
evenly throughout the year and no events are coincident 
in time with the IceTop triggers. Figure 6 shows the 
energy spectrum. The average reconstructed energy is 
{E) ~1.7 TeV with 970 events (90%) between 278 GeV 
and 13.5 TeV. Above 10 TeV, 70 events are detected. Of 
the total, 57% are reconstructed as upward-going. The 
baseline fit results are shown in Table El with the total 
uncertainties. The 1-dimensional projected distributions 
for the best ht zenith angle and PID variables are shown 
in Figures 7 and 8 . 

Figure | 6 ] shows the energy spectra of the different 
components. The CR muon events (11% of 

the total) follow the distribution estimated from the 
background-dominated data region. At the lowest en¬ 
ergies, the background contamination is track-like, and 
mostly downward-going. These events are not easily de¬ 
tectable by visual inspection and are most likely stealth 
muon events with little veto information. 

The conventional and Ve components are mainly 
horizontal, at low energies. The energy and zenith angle 
distributions are similar, so much of the Vf^jv^ separation 
power comes from the PID observable. The PID sepa¬ 
rates CC from other events which have no trace of 
a muon track : NC events and Vf. CC. The Vf. CC and 
NC events are indistinguishable. The fit finds the 
normalization at x modified Honda (645 events) 

and the Vf. normalization at I.3 ^q'3 x modified Honda 
(215 events). The flux ratio of to v^. at 1.7 TeV is 
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FIG. 10. Best fit contours for the conventional flux at 68 % 
C.L. The baseline fit with the astrophysical spectral index 7 
free is shown with solid blue line. An alternative fit with the 
index fixed at 7 = 2.0 is shown as a dotted line. 
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16.97;^ compared with the Honda prediction of 23 and 
the Bartol [l2| prediction of 14. The models use different 
assumptions about the primary cosmic-ray spectrum and 
the treatment of kaons 6 l|, l62| . 

The to ratio depends on the kaon to pion ratio in 
cosmic-ray air showers. One of the major uncertainties 
in the K : tt ratio is due to associated production via 
reactions like p -I- V — >■ A -|- 77+. A higher rate of associ¬ 
ated production leads to fewer X/g and more zzg at energies 
above I TeV . Since the Pg and iZg have different inter¬ 
action cross-sections in the ice, this will lead to a smaller 
amount in the total r'g rate, resulting in higher r'^/t'g 
ratio. Both calculations suffer from large uncertainties 
regarding kaon production at these energies. 

The statistical uncertainties on the and z/g normal¬ 
izations are estimated to be 8 . 6 % and 20 %, respectively 
as determined by running the fitter without the system¬ 
atic parameters included. The conventional normaliza¬ 
tion results are consistent with Honda predictions and 
the significance contours of the conventional normaliza¬ 
tion fit are shown in Fig. (TOl Overall, the CR muons 
and the conventional neutrinos are not correlated with 
prompt or astrophysical components. As can be seen in 
Fig. 1101 the change in conventional normalization with 
the astrophysical model is minimal. 

On the other hand, the prompt normalization is 
strongly influenced by astrophysical models. The fit 
for the prompt normalization is zero with the 68 % 
confidence upper limit at 3.0xmodified ERS. The best 
fit astrophysical flux per flavor is 3.2+7 x 10 “^® 
GeV-icm-2sr-is-i(E7l0®GeV)-'^ with 7 = 2.4l°+. 
The relationship between the fits for the prompt flux and 
astrophysical models is shown in Fig. [TTJ As the astro- 
physical spectral index softens, the shapes of the prompt 
and astrophysical components in the observable space be¬ 
come similar. In the limit of identical indices, the main 
way to separate these two components is via self-vetoing; 
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Astrophysical spectral index (y) 


FIG. 11. For a given astrophysical spectral index (x axis) 
in the upper panel, the best fit prompt flux (blue line) and 
its errors (band at 68% C.L.) from the profile likelihood scan 
are obtained. The bottom panel shows the range of allowed 
region of the index parameter from the full ht. 


down-going prompt neutrinos will be accompanied by 
muons which will cause the event to be rejected. This 
will show up as a change in the zenith angle distribution, 
with down-going events suppressed, in contrast to the 
astrophysical component, which will remain isotropic. 

The presence of very high energy events ('^1 PeV) in 
the downward region favors the astrophysical component 
over the prompt component. It should be noted that the 
presence of the cosmic-ray knee introduces a kink into 
the prompt component spectrum. As Fig. [1^ shows, 
at energies above a few hundred TeV, this kink further 
reduces the prompt component. 

Since the fit results for the conventional components 
are not influenced by the prompt or astrophysical com¬ 
ponents, we obtain the conventional spectrum inde¬ 
pendent of assumptions about the other components. A 
separate fit is performed by introducing conventional t'e 
components divided into four true energy ranges while 
keeping all of the other components unchanged. The re¬ 
sulting best-fit normalizations in each range produce the 
neutrino fluxes as shown in Fig. [12] and Table nni The 
fit finds good agreement with models of the conventional 
Vf. flux. The other components in the fit show consistent 
values when compared to the previous baseline fit. 

The relatively high conventional v^. flux normalization 
measured in the first fit can be further examined by vary¬ 
ing the relative contribution from tt and K to the con¬ 
ventional neutrino fluxes. In a third fit, we introduce an 



FIG. 12. The atmospheric flux result (shown as red filled 
triangles). Markers indicate the IceCube measurements of 
the atmospheric neutrino flux while lines show the theoreti¬ 
cal models. The black circles and the blue band come from 
the through-going upward analyses ii. The open tri¬ 
angles show the Ve measurement with the IceCube-DeepCore 
dataset [ 2 ]. The magenta band shows the modified ERS pre¬ 
diction. 

TABLE III. The results of the binned (‘second’) fit to the Ve 
flux for an E~^ spectrum, in four energy bins. 


logioEr"-logioE“- 

(E.)(GeV) 

Ej-l-,(GeV cm-^s-^sr-®) 

2.0 - 2.5 

270 

(1.0 ±0.9) X 10“® 

2.5 - 3.0 

590 

(7.6 ± 1.9) X 10“® 

3.0 - 4.0 

2.5 X 10® 

(6.4 ± 2.6) X 10“'^ 

4.0 - 5.0 

20.7 X 10® 

(3.5 ±3.3) X 10“® 


extra fit parameter (^) which modihes the K contribu¬ 
tions in Eq. |7| and in Eq. |8| simultaneously. 

(7) 

= ( 8 ) 

A value of ^ = 1 corresponds to the standard expec¬ 
tations based on the modified Honda model and a value 
of ^ > 1 corresponds to increased kaon production. As 
the conventional and flux normalizations are fixed 
to the baseline model, ^ probes the deviations from the 
model due to relative K contribution. The r'e normaliza¬ 
tion C and the kaon weight wk' are fixed at the Honda 
flux. For the part, while the change in ^ corresponds 
to a change in shape of the energy distribution, the total 
number of events is fixed to the baseline expectation 
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due to the change in On the other hand, an increase 
in the K contribution to i>e causes the number of events 
in the prediction to increase while the shape is un¬ 
changed. This is because comes mostly from K in 
these energies. The Ve flux from tt —^ decays is 
negligible, so there is little shape change in the r'e energy 
spectrum due to tt. This fit finds ^ with respect 

to the modified Honda flux. 

The central value of the K content is above stan¬ 
dard calculations, although the errors are large. Cur¬ 
rent models of cosmic ray interactions may underestimate 
the strange quark content in the air shower. Enhanced 
strangeness production has been measured in nuclear col¬ 
lisions at Relativistic Heavy Ion Collider 6^, and air 
shower experiments also measure higher muon contents 
for inclined showers compared with t he p redictions from 
existing hadronic interaction models [65l - l67j . 


IX. CONCLUSIONS 

In conclusion, we obtained a sample of 1078 cascade 
events in the analysis of one year of data from the com¬ 
pleted IceCube detector. This sample is used to measure 
the conventional atmospheric Vg flux. The analysis is de¬ 
signed so that the conventional neutrino result is largely 
unaffected by the prompt neutrino flux and/or the astro- 
physical models. The analysis extends previous measure¬ 
ments Q of the Vg flux to higher energies, and provides 
higher precision. The first analysis with only the Deep- 
Core region as a fiducial volume was optimized in obtain¬ 
ing a large number of lower energy events. Therefore, the 
improvement comes from a better event selection by ex¬ 
panding the fiducial volume for higher energy events and 
a three dimensional likelihood method including particle 
identification at higher energies. 

The conventional Vg spectrum was measured between 
0.1 TeV and 100 TeV. The measured Vg flux was 1.3to;3X 
modified Honda prediction which includes a model of 
the cosmic-ray knee and a correction to account for self- 
vetoing, whereby an atmospheric neutrino is accompa¬ 
nied by muons from the same shower, causing it to fail 
the event selection. An unfolding was used to determine 
the Vg flux in four energy bins. 

In addition to the conventional Vg spectrum measure¬ 
ments, we find that the result for the prompt component 
strongly depends on the assumed astrophysical models. 
The analysis fits the prompt flux at O.OOIq'qX modified 
ERS, together with the astrophysical flux per flavor at 
X 10“^® GeV-icm-2sr-is-i(E,,/10®GeV)-'^ with 
7 = 2.4to;2 C.L. The uniqueness of the prompt 

compared to soft astrophysical components is two-fold: a 
shape difference in energy due to the presence of cosmic 
ray knee and a shape difference in zenith angle due to 
the impact of the self veto. 

The analysis also finds a slightly higher K contribu¬ 


tion than in current models, at 1.3lg'4X modified Honda. 
The measured neutrino flux ratio v^jvg = g at the 

mean neutrino energy of 1.7 TeV, is below the prediction 
of the Honda model, but slightly above the prediction of 
the Bartol model. 

At energies above a few TeV, additional data, as would 
be provided by a multi-year analysis, would allow for a 
more precise measurement. 
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